Palmate Newts Rum Analysis 2025
Xav Harrison & Emma Remotti
2023-03-12
Data and Libraries
Libraries
library(phyloseq) #if not installed: install.packages('devtools')
library(devtools) #if not installed: install_github("microbiome/microbiome") # Install the package## Loading required package: usethis library(microbiome)## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.5.2##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>##
## Attaching package: 'microbiome'## The following object is masked from 'package:ggplot2':
##
## alpha## The following object is masked from 'package:base':
##
## transform library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union library(ggplot2)
library(cowplot)
library(RColorBrewer)
library(remotes) #if not installed: remotes::install_github("jfq3/ggordiplots")##
## Attaching package: 'remotes'## The following objects are masked from 'package:devtools':
##
## dev_package_deps, install_bioc, install_bitbucket, install_cran,
## install_deps, install_dev, install_git, install_github,
## install_gitlab, install_local, install_svn, install_url,
## install_version, update_packages library(ggordiplots)## Loading required package: vegan## Loading required package: permute##
## Attaching package: 'permute'## The following object is masked from 'package:devtools':
##
## check##
## Attaching package: 'vegan'## The following object is masked from 'package:microbiome':
##
## diversity## Loading required package: glue library(vegan)
library(ALDEx2) #to install: BiocManager::install("ALDEx2")## Loading required package: zCompositions## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: truncnorm## Loading required package: survival## Loading required package: lattice## Loading required package: latticeExtra##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layer library(decontam)
library(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.23.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:survival':
##
## kidney## The following object is masked from 'package:phyloseq':
##
## nsamples## The following object is masked from 'package:stats':
##
## ar library(MoMAColors)
library(ggrepel)
library(ggvenn)
library(gllvm)## Loading required package: TMB##
## Attaching package: 'gllvm'## The following object is masked from 'package:vegan':
##
## ordiplot## The following object is masked from 'package:stats':
##
## simulate library(microbiomeutilities) #devtools::install_github("microsud/microbiomeutilities")## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'##
## Attaching package: 'microbiomeutilities'## The following object is masked from 'package:microbiome':
##
## add_refseq library(ggvegan)Global Plotting Options
Quick access options for output graphics to make legends and axis labels larger / more legible
### Site Colours for Visit One From Paired Palette
sitecols<-brewer.pal(8,"Paired")[c(1,3,5,7)]Data
####### READ IN FILES
#Taxonomy
tax1<-readRDS('Newt 2021 Taxonomy.rds')
#rownames(tax1)<-paste0("SV",seq(nrow(tax1)))
#Sample Metadata
sample1<-readRDS('Newt 2021 SampleData.rds')
head(sample1)## Sample
## 1 DP-1
## 2 DP-10
## 3 DP-11
## 4 DP-12
## 5 DP-13
## 6 DP-14 colnames(sample1)[1]<-"swab"
sample1$visit<-ifelse(grepl("C$",sample1$swab),2,1)
sample1$swab_time<-with(sample1,paste(gsub("C$","",swab),visit,sep="_"))
write.csv(sample1,'Newt Sample Data.csv')
## Additional Metadata
meta1<-read.csv('Newt Metadata Full Feb17 R Input.csv')
head(meta1)## visit airtemp watertemp ph sex svl infection swab
## 1 1 12.3 13.1 5.11 1 4.0 0 DP-1
## 2 1 12.3 13.1 5.11 1 3.9 0 DP-2
## 3 1 12.3 13.1 5.11 1 3.5 0 DP-3
## 4 1 12.3 13.1 5.11 1 3.6 0 DP-4
## 5 1 12.3 13.1 5.11 F 3.9 0 DP-5
## 6 1 12.3 13.1 5.11 F 4.0 0 DP-6 meta1$swab_time<-with(meta1,paste(swab,visit,sep="_"))
#Copy Across Metadata
sample2<-sample1 %>% dplyr::select(swab_time) %>% left_join(meta1,"swab_time")
#head(sample2)
sample2$site<- sapply(strsplit(sample2$swab,"-"),"[",1)
rownames(sample2)<-sample1[,1]
##Sample Type
sample2$sampletype<-"swab"
sample2$sampletype[grep("MC",rownames(sample2))]<-"mock"
sample2$sampletype[grep("NEG",rownames(sample2))]<-"negative_control"
#Sequence Abundance
otu1<-readRDS('Newt 2021 SeqTab.rds')
#head(otu1)
#Infection Data
infection<-read.csv('newt infection codes.csv',header=T)
head(infection)## swab lesion timecode
## 1 DP-1 <NA> 1
## 2 DP-2 <NA> 1
## 3 DP-3 <NA> 1
## 4 DP-4 <NA> 1
## 5 DP-5 <NA> 1
## 6 DP-6 <NA> 1 infection$swab_time<-with(infection,paste(swab,timecode,sep="_"))
sample2$infectioncode<-infection$lesion[match(sample2$swab_time,infection$swab_time)]
head(sample2)## swab_time visit airtemp watertemp ph sex svl infection swab site
## DP-1 DP-1_1 1 12.3 13.1 5.11 1 4.0 0 DP-1 DP
## DP-10 DP-10_1 1 12.3 13.1 5.11 1 3.5 0 DP-10 DP
## DP-11 DP-11_1 1 12.3 13.1 5.11 1 3.5 0 DP-11 DP
## DP-12 DP-12_1 1 12.3 13.1 5.11 1 3.3 0 DP-12 DP
## DP-13 DP-13_1 1 12.3 13.1 5.11 1 3.6 0 DP-13 DP
## DP-14 DP-14_1 1 12.3 13.1 5.11 1 3.9 0 DP-14 DP
## sampletype infectioncode
## DP-1 swab <NA>
## DP-10 swab <NA>
## DP-11 swab <NA>
## DP-12 swab <NA>
## DP-13 swab <NA>
## DP-14 swab <NA> #Re-Name Rows with COrrect Sample IDs
rownames(otu1)<-sample1[,1]Phyloseq Build
newt <- phyloseq(otu_table(otu1,taxa_are_rows = FALSE), # create new phyloseq otu table, call it "newt"
sample_data(sample2),
tax_table((tax1)))
newt ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 13903 taxa and 152 samples ]
## sample_data() Sample Data: [ 152 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 13903 taxa by 7 taxonomic ranks ]Decontamination
##### REMOVE MOCK SAMPLES
newt.nomock<-prune_samples(sample_data(newt)$sampletype %in% c("swab","negative_control"),newt)
newt.nomock## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 13903 taxa and 150 samples ]
## sample_data() Sample Data: [ 150 samples by 12 sample variables ]
## tax_table() Taxonomy Table: [ 13903 taxa by 7 taxonomic ranks ] #colnames(tax_table(newt.nomock))<-c("Kingdom","Phylum","Class","Order","Family","Genus","Species")
###### REMOVE NEGATIVE CONTAMINATION
sample_data(newt.nomock)$is.neg <- ifelse(grepl("negative_control",sample_data(newt.nomock)$sampletype),TRUE,FALSE)
contamdf.prev <- isContaminant(newt.nomock, method="prevalence", neg="is.neg",threshold=0.5)
table(contamdf.prev$contaminant)##
## FALSE TRUE
## 13870 33 # 33 contaminants identified
###### Plot of Contaminant Frequency
# Make phyloseq object of presence-absence in negative controls and true samples
ps.pa <- transform_sample_counts(newt.nomock, function(abund) 1*(abund>0))
ps.pa.neg <- prune_samples(sample_data(ps.pa)$is.neg == TRUE, ps.pa)
ps.pa.pos <- prune_samples(sample_data(ps.pa)$is.neg == FALSE, ps.pa)
# Make data.frame of prevalence in positive and negative samples
df.pa <- data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg),
contaminant=contamdf.prev$contaminant)
ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() +
xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")
###### Clean Contaminants Out
newt.clean<-prune_taxa(contamdf.prev$contaminant==FALSE,newt.nomock)
ntaxa(newt.nomock) - ntaxa(newt.clean)## [1] 33Further Cleaning
##### REMOVE NEGATIVES SAMPLES
newt.clean<-prune_samples(sample_data(newt.clean)$sampletype=="swab",newt.clean)
newt.clean## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 13870 taxa and 148 samples ]
## sample_data() Sample Data: [ 148 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 13870 taxa by 7 taxonomic ranks ] #Prune Taxa With No Phylum Assignment
ps_prune<-prune_taxa(as.vector(!is.na(tax_table(newt.clean)[,2])),newt.clean)
ntaxa(newt.clean)-ntaxa(ps_prune) #232 lost## [1] 232 #Prune Chloroplasts
ps_prune_nochloro<-prune_taxa(as.vector(tax_table(ps_prune)[,3]!="Chloroplast"),ps_prune)
ntaxa(ps_prune)-ntaxa(ps_prune_nochloro) #431 lost## [1] 431 #Remove Mitochondria
ps_prune_nochloro_nomito<-prune_taxa(as.vector(tax_table(ps_prune_nochloro)[,5]!="Mitochondria"),ps_prune_nochloro)
#Remove Archaea
ps_prune_nochloro_noarchaea<-prune_taxa(as.vector(tax_table(ps_prune_nochloro_nomito)[,1]!="Archaea"),ps_prune_nochloro_nomito)
ntaxa(ps_prune_nochloro_noarchaea)-ntaxa(ps_prune_nochloro_nomito) #93 lost## [1] -93Filter Based on Prevalence and Abundance
####### Read Threshold and Sample Prevalence Threshold
ps_clean_filter = filter_taxa(ps_prune_nochloro_noarchaea, function(x) sum(x > 20) > (0.02*length(x)), TRUE)Mock Community
##Inspect Mock After Removing ASVs filtered at previous step
ps_mock<-prune_samples(sample_data(newt)$sampletype %in% c("mock"),newt)
ps_mock<-filter_taxa(ps_mock, function(x) (sum(x) > 0),TRUE)
data.frame(reads=taxa_sums(ps_mock),genus=unname(tax_table(ps_mock))[,6])## reads
## TACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGCTTGATAAGCCGGTTGTGAAAGCCCCGGGCTCAACCTGGGAACGGCATCCGGAACTGTCAAGCTAGAGTGCAGGAGAGGAAGGTAGAATTCCCGGTGTAGCGGTGAAATGCGTAGAGATCGGGAGGAATACCAGTGGCGAAGGCGGCCTTCTGGACTGACACTGACACTGAGGTGCGAAAGCGTGGGTAGCAAACAGG 11
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG 1322
## TACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG 21387
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCAAAACTGTCGAGCTAGAGTATGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG 13464
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAATTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG 15499
## TACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATATTTAAGTCAGGGGTGAAATCCCAGAGCTCAACTCTGGAACTGCCTTTGATACTGGGTATCTTGAGTATGGAAGAGGTAAGTGGAATTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAGGAACACCAGTGGCGAAGGCGGCTTACTGGTCCATAACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG 6853
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAGTTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG 2507
## TACGTAGGGTCCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCTGGAGGCTCAACCTCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGAATGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGTTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACAGG 2248
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG 1676
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTGAGGCGCGAAAGCGTGGGGAGCAAACAGG 1790
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG 1364
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCTTAGCTTAACTAAGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG 1085
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTGATGTGCGAAAGCGTGGGGATCAAACAGG 893
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTTGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG 932
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCATGCAGGTGGTTTGTTAAGTCAGATGTGAAAGCCCTGGGCTCAACCTAGGAATCGCATTTGAAACTGACAAGCTAGAGTACTGTAGAGGGGGGTAGAATTTCAGGTGTAGCGGTGAAATGCGTAGAGATCTGAAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAGATACTGACACTCAGATGCGAAAGCGTGGGGAGCAAACAGG 857
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCATTTGAAACTGGCAAGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG 787
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG 703
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTGATTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGTCAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG 654
## TACGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTTGATTAAGTGAGATGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGTCAACTAGAGTACTTTAGGGAGGGGTAGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAATACCGAAGGCGAAGGCAGCCCCTTGGGAATGTACTGACGCTCATGTGCGAAAGCGTGGGGAGCAAACAGG 438
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTGGCGAGCTAGAGTAGGGCAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGGCTCATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG 323
## TACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAGCCCGGGGCTTAACCCCGGGTGTGCAGTGGGTACGGGCAGACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCTGTTACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACAGG 200
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGCTTTGTAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACTGCCATTGATACTGCAGAGCTTGAATAATATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTATTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG 278
## TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTGATAAGCCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCGTTTGGAACTGTCAGACTAGAGTGCGTCAGAGGGGGGTGGAATTCCGCGTGTAGCAGTGAAATGCGTAGAGATGCGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGATGACACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACAGG 76
## TACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCCCGAGGCTCAACCTCGGGTCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGATCTCTGGGCCGCTACTGACGCTGAGGAGCGAAAGGGTGGGGAGCAAACAGG 54
## CGATGAACATGATTAGCGGTACGACTTTGGCCACCGTGGTCAGTTGATTGATCAGCGCCGCCTCCTTGATCCCGCGCATCACCAAAAAATGCACGGCCCACAGCAGCACCGACGCGCAGCCGATGGCGATCGGCGTATTGCCTTGGCCAAACACCGGAAAGAAGTAGCCGAGGGTGCTGAACAGCAGGACGAAATAACCGACGTTGCCCAGCCACGCACTGATCCAGTAACCCCAGGCTG 6
## ta1
## TACGGAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGCTTGATAAGCCGGTTGTGAAAGCCCCGGGCTCAACCTGGGAACGGCATCCGGAACTGTCAAGCTAGAGTGCAGGAGAGGAAGGTAGAATTCCCGGTGTAGCGGTGAAATGCGTAGAGATCGGGAGGAATACCAGTGGCGAAGGCGGCCTTCTGGACTGACACTGACACTGAGGTGCGAAAGCGTGGGTAGCAAACAGG Halomonas
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG Klebsiella
## TACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG Sphingomonas
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCAAAACTGTCGAGCTAGAGTATGGTAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGACTGATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG Pseudomonas
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAATTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG Flavobacterium
## TACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATATTTAAGTCAGGGGTGAAATCCCAGAGCTCAACTCTGGAACTGCCTTTGATACTGGGTATCTTGAGTATGGAAGAGGTAAGTGGAATTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAGGAACACCAGTGGCGAAGGCGGCTTACTGGTCCATAACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGG Neorhizobium
## TACGGAGGATCCAAGCGTTATCCGGAATCATTGGGTTTAAAGGGTCCGTAGGCGGTTTAATAAGTCAGTGGTGAAAGCCCATCGCTCAACGGTGGAACGGCCATTGATACTGTTAGACTTGAATTATTAGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGAGATTACATGGAATACCAATTGCGAAGGCAGGTTACTACTAGTTGATTGACGCTGATGGACGAAAGCGTGGGTAGCGAACAGG Flavobacterium
## TACGTAGGGTCCGAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCTGGAGGCTCAACCTCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGAATGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGTTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACAGG Agromyces
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG Citrobacter
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGTGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAGACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGTAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACACTGAGGCGCGAAAGCGTGGGGAGCAAACAGG Bacillus
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG Myroides
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGTCAGATAAGTCAGTGGTGAAATTTCTTAGCTTAACTAAGACACGGCCATTGATACTGTTTGACTTGAATAGTATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTACTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG Myroides
## TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTGATGTGCGAAAGCGTGGGGATCAAACAGG Staphylococcus
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTTGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG Salmonella
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCATGCAGGTGGTTTGTTAAGTCAGATGTGAAAGCCCTGGGCTCAACCTAGGAATCGCATTTGAAACTGACAAGCTAGAGTACTGTAGAGGGGGGTAGAATTTCAGGTGTAGCGGTGAAATGCGTAGAGATCTGAAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAGATACTGACACTCAGATGCGAAAGCGTGGGGAGCAAACAGG Vibrio
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCATTTGAAACTGGCAAGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG Plesiomonas
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTTGTTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGCAAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG Escherichia-Shigella
## TACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTTGATTAAGTCAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCTGATACTGGTCAGCTTGAGTCTCGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACGAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGG <NA>
## TACGGAGGGTGCGAGCGTTAATCGGAATAACTGGGCGTAAAGGGCACGCAGGCGGTTGATTAAGTGAGATGTGAAAGCCCCGGGCTTAACCTGGGAATTGCATTTCATACTGGTCAACTAGAGTACTTTAGGGAGGGGTAGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAATACCGAAGGCGAAGGCAGCCCCTTGGGAATGTACTGACGCTCATGTGCGAAAGCGTGGGGAGCAAACAGG Actinobacillus
## TACAGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTTCGTTAAGTTGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATCCAAAACTGGCGAGCTAGAGTAGGGCAGAGGGTGGTGGAATTTCCTGTGTAGCGGTGAAATGCGTAGATATAGGAAGGAACACCAGTGGCGAAGGCGACCACCTGGGCTCATACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGG Pseudomonas
## TACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAGCCCGGGGCTTAACCCCGGGTGTGCAGTGGGTACGGGCAGACTAGAGTGCAGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCTGTTACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACAGG Kocuria
## TACGGAGGATCCGAGCGTTATCCGGAATTATTGGGTTTAAAGGGTTCGTAGGCGGCTTTGTAAGTCAGTGGTGAAATTTCCTAGCTTAACTAGGACACTGCCATTGATACTGCAGAGCTTGAATAATATGGAAGTAACTAGAATATGTAGTGTAGCGGTGAAATGCTTAGATATTACATGGAATACCAATTGCGAAGGCAGGTTACTACGTATTTATTGACGCTGATGAACGAAAGCGTGGGGAGCGAACAGG Myroides
## TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTGATAAGCCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCGTTTGGAACTGTCAGACTAGAGTGCGTCAGAGGGGGGTGGAATTCCGCGTGTAGCAGTGAAATGCGTAGAGATGCGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGATGACACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACAGG Vogesella
## TACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAATCCCGAGGCTCAACCTCGGGTCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGATCTCTGGGCCGCTACTGACGCTGAGGAGCGAAAGGGTGGGGAGCAAACAGG Plantibacter
## CGATGAACATGATTAGCGGTACGACTTTGGCCACCGTGGTCAGTTGATTGATCAGCGCCGCCTCCTTGATCCCGCGCATCACCAAAAAATGCACGGCCCACAGCAGCACCGACGCGCAGCCGATGGCGATCGGCGTATTGCCTTGGCCAAACACCGGAAAGAAGTAGCCGAGGGTGCTGAACAGCAGGACGAAATAACCGACGTTGCCCAGCCACGCACTGATCCAGTAACCCCAGGCTG <NA>##Library Sizes After Cleaning
#Summary
range(sample_sums(ps_clean_filter)) ## [1] 1278 69644 mean(sample_sums(ps_clean_filter)) ## [1] 22264.95#Table
sample_sums(ps_clean_filter)## DP-1 DP-10 DP-11 DP-12 DP-13 DP-14 DP-15 DP-16 DP-17 DP-18
## 28480 38000 33723 36636 56889 42389 32158 34777 27873 28435
## DP-19 DP-1C DP-2 DP-20 DP-21 DP-22 DP-23 DP-24 DP-25 DP-26
## 32201 14498 43533 32864 39365 36527 31885 42616 52852 43414
## DP-27 DP-2C DP-3 DP-3C DP-4 DP-4C DP-5 DP-5C DP-6 DP-7
## 27425 25949 35236 17054 69544 20893 53744 15085 38979 55027
## DP-8 DP-8C DP-9 R10-1 R10-10 R10-11 R10-12 R10-13 R10-14 R10-15
## 69644 14214 24015 11275 8525 16545 14457 9632 17585 14677
## R10-16 R10-17 R10-18 R10-19 R10-1C R10-2 R10-20 R10-2C R10-3 R10-3C
## 19851 12916 11771 11105 15820 15172 15998 28945 9717 28280
## R10-4 R10-5 R10-6 R10-6C R10-7 R10-7C R10-8 R10-8C R4-10 R4-11
## 9688 10696 13754 9088 9817 18480 13746 19026 7882 26514
## R4-12 R4-12C R4-13 R4-13C R4-14 R4-14C R4-18C R4-1C R4-2 R4-21
## 10092 21847 20223 17351 15567 23889 9373 11501 16158 27400
## R4-21C R4-22 R4-22C R4-23 R4-23C R4-24 R4-25 R4-26 R4-27 R4-28
## 21423 7133 16220 19078 18841 9133 12559 13181 6693 19605
## R4-2C R4-3 R4-30 R4-31 R4-32 R4-33 R4-34 R4-35 R4-36 R4-38
## 8476 28765 9835 7668 1278 20874 13017 7076 29156 24914
## R4-39 R4-3C R4-4 R4-40 R4-41 R4-42 R4-43 R4-44 R4-45 R4-45C
## 22101 13156 51135 14121 11672 10043 18843 31806 10400 27929
## R4-46 R4-47 R4-48 R4-49 R4-4C R4-50 R4-5C R4-6 R4-6C R4-7
## 22375 45324 5989 24170 18693 24692 30438 27978 16159 21036
## R4-8 R4-9 R44-1 Rnew-1 Rnew-10 Rnew-11 Rnew-12 Rnew-13 Rnew-14 Rnew-15
## 15286 16207 11206 17843 18919 19913 28563 21616 23889 36707
## Rnew-2 Rnew-3 Rnew-37 Rnew-38 Rnew-39 Rnew-4 Rnew-40 Rnew-41 Rnew-42 Rnew-43
## 28399 41446 23936 23781 8069 27542 8690 21795 25391 13301
## Rnew-44 Rnew-45 Rnew-46 Rnew-47 Rnew-48 Rnew-49 Rnew-5 Rnew-50 Rnew-51 Rnew-52
## 21870 14670 5943 8300 7312 12217 49985 11141 5390 8642
## Rnew-53 Rnew-54 Rnew-55 Rnew-56 Rnew-6 Rnew-7 Rnew-8 Rnew-9
## 7532 8148 20280 19401 38318 18771 52717 38805Prevalence Models
Probability of Observing C and D Infection
sample2$C_or_D<-ifelse(sample2$infectioncode %in% c("C","D"),1,0)
#table(sample2$infectioncode)
with(sample2,table(infectioncode,site))## site
## infectioncode DP R10 R4 Rnew
## A 0 0 11 10
## B 0 0 14 6
## C 0 0 9 4
## D 0 0 2 0prev_mod1<-brm(C_or_D ~ site,data = sample2,family=bernoulli()) ## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.12 seconds (Warm-up)
## Chain 1: 0.126 seconds (Sampling)
## Chain 1: 0.246 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.124 seconds (Warm-up)
## Chain 2: 0.102 seconds (Sampling)
## Chain 2: 0.226 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.109 seconds (Warm-up)
## Chain 3: 0.146 seconds (Sampling)
## Chain 3: 0.255 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.094 seconds (Warm-up)
## Chain 4: 0.109 seconds (Sampling)
## Chain 4: 0.203 seconds (Total)
## Chain 4:## Warning: There were 2 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess summary(prev_mod1)## Family: bernoulli
## Links: mu = logit
## Formula: C_or_D ~ site
## Data: sample2 (Number of observations: 147)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -21.33 34.53 -92.28 -3.67 1.02 145 133
## siteR10 -10.76 67.42 -128.74 43.14 1.02 232 126
## siteR4 19.95 34.53 2.27 91.20 1.02 144 134
## siteRnew 19.20 34.53 1.46 90.16 1.02 146 133
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).ALPHA DIVERSITY
Alpha Calculations First Visit
##Rarefy
newt1<-prune_samples(sample_data(ps_clean_filter)$visit==1 & sample_data(ps_clean_filter)$site!="R44",ps_clean_filter)
newt1 #121 samples after removing R44## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 893 taxa and 121 samples ]
## sample_data() Sample Data: [ 121 samples by 13 sample variables ]
## tax_table() Taxonomy Table: [ 893 taxa by 7 taxonomic ranks ]#Rarefy data 5000 reads
newt1r<-rarefy_even_depth(newt1,rngseed = 01042020,sample.size = 5000) ## `set.seed(1042020)` was used to initialize repeatable random subsampling.## Please record this for your records so others can reproduce.## Try `set.seed(1042020); .Random.seed` for the full vector## ...## 1 samples removedbecause they contained fewer reads than `sample.size`.## Up to first five removed samples are:## R4-32## ...## 6OTUs were removed because they are no longer
## present in any sample after random subsampling## ... #range(sample_sums(newt1r))
##Estimate alpha diversity metrics
newt1_rich<-estimate_richness(newt1r,measures=c("Observed"))
newt1_rich$swab<-gsub("\\.","-",rownames(newt1_rich))
newt1_rich$swab_time<-with(newt1_rich,paste(swab, "1",sep="_"))
##Join
newt1_rich_meta<-left_join(newt1_rich,sample2,"swab_time")
head(newt1_rich_meta)## Observed swab.x swab_time visit airtemp watertemp ph sex svl infection
## 1 81 DP-1 DP-1_1 1 12.3 13.1 5.11 1 4.0 0
## 2 112 DP-10 DP-10_1 1 12.3 13.1 5.11 1 3.5 0
## 3 143 DP-11 DP-11_1 1 12.3 13.1 5.11 1 3.5 0
## 4 55 DP-12 DP-12_1 1 12.3 13.1 5.11 1 3.3 0
## 5 111 DP-13 DP-13_1 1 12.3 13.1 5.11 1 3.6 0
## 6 88 DP-14 DP-14_1 1 12.3 13.1 5.11 1 3.9 0
## swab.y site sampletype infectioncode C_or_D
## 1 DP-1 DP swab <NA> 0
## 2 DP-10 DP swab <NA> 0
## 3 DP-11 DP swab <NA> 0
## 4 DP-12 DP swab <NA> 0
## 5 DP-13 DP swab <NA> 0
## 6 DP-14 DP swab <NA> 0##Standardise and Factorise
newt1_rich_meta$z_svl<- with(newt1_rich_meta,as.numeric(scale(svl))) # mean centered svl (snout-ventral length)
newt1_rich_meta$watertemp<-factor(newt1_rich_meta$watertemp)
newt1_rich_meta$infection<-factor(newt1_rich_meta$infection)
#newt1_rich_meta$ph<-factor(newt1_rich_meta$ph)
newt1_rich_meta$site<-factor(newt1_rich_meta$site)Alpha Models
Site, Sex, SVL
#Maximal Model
rich_mod1<-brm(Observed ~ site + sex + z_svl, data = newt1_rich_meta, family = negbinomial())## Compiling Stan program...## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.07 seconds (Warm-up)
## Chain 1: 0.074 seconds (Sampling)
## Chain 1: 0.144 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.071 seconds (Warm-up)
## Chain 2: 0.068 seconds (Sampling)
## Chain 2: 0.139 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.071 seconds (Warm-up)
## Chain 3: 0.068 seconds (Sampling)
## Chain 3: 0.139 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.068 seconds (Warm-up)
## Chain 4: 0.067 seconds (Sampling)
## Chain 4: 0.135 seconds (Total)
## Chain 4: summary(rich_mod1)## Family: negbinomial
## Links: mu = log
## Formula: Observed ~ site + sex + z_svl
## Data: newt1_rich_meta (Number of observations: 120)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.58 0.07 4.44 4.72 1.00 3517 2815
## siteR10 0.70 0.14 0.44 0.97 1.00 2573 2392
## siteR4 -0.55 0.09 -0.72 -0.38 1.00 3011 2881
## siteRnew -0.74 0.09 -0.92 -0.55 1.00 2801 2769
## sexF -0.02 0.09 -0.20 0.16 1.00 2992 2249
## sexM -0.19 0.18 -0.55 0.15 1.00 2673 2437
## z_svl -0.02 0.04 -0.10 0.07 1.00 3060 2836
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 9.53 1.44 6.96 12.57 1.00 3495 3207
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).#Diagnostics
pp_check(rich_mod1)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
bayes_R2(rich_mod1)## Estimate Est.Error Q2.5 Q97.5
## R2 0.7427242 0.03864956 0.6462914 0.7939698#Plot
conditional_effects(rich_mod1)## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
#Null Model
rich_mod_null<-brm(Observed ~ 1, data = newt1_rich_meta, family = negbinomial())## Compiling Stan program...
## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.045 seconds (Warm-up)
## Chain 1: 0.041 seconds (Sampling)
## Chain 1: 0.086 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.044 seconds (Warm-up)
## Chain 2: 0.047 seconds (Sampling)
## Chain 2: 0.091 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 9e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.045 seconds (Warm-up)
## Chain 3: 0.051 seconds (Sampling)
## Chain 3: 0.096 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 9e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.044 seconds (Warm-up)
## Chain 4: 0.047 seconds (Sampling)
## Chain 4: 0.091 seconds (Total)
## Chain 4:#Add ICs
rich_mod1<-add_criterion(rich_mod1,"loo")
rich_mod_null<-add_criterion(rich_mod_null,"loo")
#Compare
loo_compare(rich_mod1,rich_mod_null)## elpd_diff se_diff
## rich_mod1 0.0 0.0
## rich_mod_null -65.3 9.4#Extract Estimates
rich_predict<- conditional_effects(rich_mod1)[[1]]pH
#Numeric
#newt1_rich_meta$ph<-as.numeric(as.character(newt1_rich_meta$ph))
#Model
ph_model1<-brm(Observed ~ ph + I(ph^2) + infection,data=newt1_rich_meta,family=negbinomial())## Compiling Stan program...## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.739 seconds (Warm-up)
## Chain 1: 0.695 seconds (Sampling)
## Chain 1: 1.434 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 7e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.659 seconds (Warm-up)
## Chain 2: 0.642 seconds (Sampling)
## Chain 2: 1.301 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.669 seconds (Warm-up)
## Chain 3: 0.65 seconds (Sampling)
## Chain 3: 1.319 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.689 seconds (Warm-up)
## Chain 4: 0.791 seconds (Sampling)
## Chain 4: 1.48 seconds (Total)
## Chain 4: summary(ph_model1)## Family: negbinomial
## Links: mu = log
## Formula: Observed ~ ph + I(ph^2) + infection
## Data: newt1_rich_meta (Number of observations: 120)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 41.14 4.08 33.15 49.28 1.00 1724 1834
## ph -12.89 1.47 -15.82 -9.98 1.00 1732 1883
## IphE2 1.12 0.13 0.86 1.38 1.00 1746 1841
## infection1 -0.24 0.09 -0.41 -0.08 1.00 2491 2341
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 8.72 1.31 6.36 11.46 1.00 3065 1986
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1). #Diagnostics
pp_check(ph_model1)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
bayes_R2(ph_model1)## Estimate Est.Error Q2.5 Q97.5
## R2 0.7385453 0.03498745 0.6452589 0.7784395 #Plots
conditional_effects(ph_model1)## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
#Model Selection
ph_model_null<-brm(Observed ~ 1,data=newt1_rich_meta,family=negbinomial())## Compiling Stan program...
## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.043 seconds (Warm-up)
## Chain 1: 0.046 seconds (Sampling)
## Chain 1: 0.089 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.044 seconds (Warm-up)
## Chain 2: 0.038 seconds (Sampling)
## Chain 2: 0.082 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 8e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.045 seconds (Warm-up)
## Chain 3: 0.043 seconds (Sampling)
## Chain 3: 0.088 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 8e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.047 seconds (Warm-up)
## Chain 4: 0.048 seconds (Sampling)
## Chain 4: 0.095 seconds (Total)
## Chain 4: ph_model1<-add_criterion(ph_model1,"loo")
ph_model_null<-add_criterion(ph_model_null,"loo")
loo_compare(ph_model_null,ph_model1)## elpd_diff se_diff
## ph_model1 0.0 0.0
## ph_model_null -62.2 8.7 #Extract
ph_predict<-conditional_effects(ph_model1)[[1]]Infection Class
#Subset Data
newt1_rich_infection<- newt1_rich_meta %>% mutate(infection=factor(infection)) %>% filter(site %in% c("R4","Rnew"))
#Model
infect_mod1<-brm(Observed ~ infection*site,family=negbinomial(), data=newt1_rich_infection)## Compiling Stan program...## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.044 seconds (Warm-up)
## Chain 1: 0.041 seconds (Sampling)
## Chain 1: 0.085 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.044 seconds (Warm-up)
## Chain 2: 0.034 seconds (Sampling)
## Chain 2: 0.078 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.788 seconds (Warm-up)
## Chain 3: 3.151 seconds (Sampling)
## Chain 3: 5.939 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 4e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.044 seconds (Warm-up)
## Chain 4: 0.037 seconds (Sampling)
## Chain 4: 0.081 seconds (Total)
## Chain 4:## Warning: There were 133 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 867 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is 1.59, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess summary(infect_mod1)## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.## Warning: There were 133 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: negbinomial
## Links: mu = log
## Formula: Observed ~ infection * site
## Data: newt1_rich_infection (Number of observations: 74)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -96.94 175.23 -400.40 4.42 1.53 7 25
## infection1 -10.78 18.16 -42.23 -0.06 1.56 7 14
## siteRnew 27.06 47.63 -0.71 109.55 1.55 7 11
## infection1:siteRnew 20.19 34.34 0.02 79.66 1.59 7 13
##
## Further Distributional Parameters:
## Estimate
## shape 82860327994598638985802098729997548010686563506201240169762507419768653154650947584.00
## Est.Error
## shape 143536241170127931064739043427420160456027114437337950929267711863194371876355833856.00
## l-95% CI
## shape 5.58
## u-95% CI
## shape 331441311978395688260968415246861568057379908849293253439115776383118818970234257408.00
## Rhat Bulk_ESS Tail_ESS
## shape 1.52 7 NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1). pp_check(infect_mod1)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
bayes_R2(infect_mod1)## Estimate Est.Error Q2.5 Q97.5
## R2 0.1111109 0.0894021 3.163101e-224 0.2911242 conditional_effects(infect_mod1)## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
## Null Model
infect_mod_null<-brm(Observed ~ 1,family=negbinomial(), data=newt1_rich_infection)## Compiling Stan program...
## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 1.8e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 1: 0.031 seconds (Sampling)
## Chain 1: 0.064 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.032 seconds (Warm-up)
## Chain 2: 0.035 seconds (Sampling)
## Chain 2: 0.067 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 6e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.048 seconds (Warm-up)
## Chain 3: 0.031 seconds (Sampling)
## Chain 3: 0.079 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 6e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.033 seconds (Warm-up)
## Chain 4: 0.027 seconds (Sampling)
## Chain 4: 0.06 seconds (Total)
## Chain 4: #Add ICs
infect_mod1<-add_criterion(infect_mod1,"loo")
infect_mod_null<-add_criterion(infect_mod_null,"loo")
#Select
print(loo_compare(infect_mod1,infect_mod_null),simplify=F)## elpd_diff se_diff elpd_loo se_elpd_loo p_loo
## infect_mod_null 0.0 0.0 -325.5 7.8 2.1
## infect_mod1 -1370526.2 85788.1 -1370851.8 85794.0 1370512.0
## se_p_loo looic se_looic
## infect_mod_null 0.6 651.0 15.6
## infect_mod1 85789.7 2741703.5 171588.0#Model Predictions
infect_pred_data<-data.frame(site=rep(c("R4","Rnew"),2),infection=rep(c(0,1),each=2))
infect_predict<-conditional_effects(infect_mod1,"site:infection")[[1]]Infection Severity
#combine infection categories c/d
newt1_rich_meta$infection_collapse<-as.factor(newt1_rich_meta$infectioncode)
table(newt1_rich_meta$infection_collapse)##
## A B C D
## 19 15 11 2 levels(newt1_rich_meta$infection_collapse)<-c("A","B","C","C")
#Filter
newt1_rich_meta_filter<- filter(newt1_rich_meta,!is.na(infection_collapse))
nrow(newt1_rich_meta_filter) #47## [1] 47#Infection Model
severity_model1<-brm(Observed ~ infection_collapse,data= newt1_rich_meta_filter,family=negbinomial())## Compiling Stan program...## Start sampling##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.024 seconds (Warm-up)
## Chain 1: 0.023 seconds (Sampling)
## Chain 1: 0.047 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 4e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.022 seconds (Warm-up)
## Chain 2: 0.021 seconds (Sampling)
## Chain 2: 0.043 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 4e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.023 seconds (Warm-up)
## Chain 3: 0.022 seconds (Sampling)
## Chain 3: 0.045 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.022 seconds (Warm-up)
## Chain 4: 0.021 seconds (Sampling)
## Chain 4: 0.043 seconds (Total)
## Chain 4: summary(severity_model1)## Family: negbinomial
## Links: mu = log
## Formula: Observed ~ infection_collapse
## Data: newt1_rich_meta_filter (Number of observations: 47)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 3.96 0.09 3.78 4.14 1.00 3998 2518
## infection_collapseB -0.03 0.14 -0.29 0.25 1.00 3815 2831
## infection_collapseC -0.23 0.14 -0.52 0.06 1.00 4226 2941
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 7.68 1.83 4.72 11.82 1.00 3306 2677
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1). #Diagnostics
pp_check(severity_model1)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
bayes_R2(severity_model1)## Estimate Est.Error Q2.5 Q97.5
## R2 0.07962835 0.05629761 0.003244018 0.2120769 #Plot
conditional_effects(severity_model1)## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"## Ignoring unknown labels:
## • fill : "NA"
## • colour : "NA"
Alpha Plots
Richness by Site
##Richness By Site
rich_site1 <- ggplot() + geom_jitter(data=newt1_rich_meta,aes(x=site,y=Observed, fill = site),color="white",shape=21,size=5,width=0.2) + scale_fill_brewer(palette = "Paired") + theme_bw(base_size=15) + theme(text = element_text(size = 18)) #+ guides(fill="none")
rich_site2<-rich_site1 + geom_errorbar(data=rich_predict,aes(x=site,ymin=lower__,ymax=upper__),width=0.1,lwd=1.2) + geom_point(data=rich_predict,aes(x=site,y=estimate__),shape=23,size=5,fill="white") + labs(x="Site",y="Alpha Diversity \n (Richness)") + scale_fill_manual(values=sitecols)## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.rich_site3<- rich_site2 +
theme(legend.key = element_rect(fill = "white"), legend.text = element_text(color = "white"), legend.title = element_text(color = "white")) + guides(fill = guide_legend(override.aes = list(fill = NA)))
rich_site3
ggsave2('Richness by Site Dec23.pdf',rich_site3,height=15,width=16,units="cm")Shannon Split by Infection Status
rich_infection<- ggplot() + geom_point(data=newt1_rich_infection,aes(x=site,y=Observed,fill=site,shape=infection),position = position_jitterdodge(),size=5,color="white") + theme_bw(base_size=15) + theme(legend.position = "right") + scale_shape_manual(name="Amphibiothecum \n Infection",values=c(21,24),labels=c("No Visible \n Infection","Infected"))
rich_infection2<- rich_infection + labs(y="Alpha Diversity \n (Richness)",x="Site") #+ theme(legend.text = element_text(size=12),legend.title = element_text(size=12))
rich_infection3<- rich_infection2 + geom_errorbar(data=infect_predict,aes(x=site,ymin=lower__,ymax=upper__,group=infection),position = position_dodge(width=0.7),width=0.1,lwd=1.2) + geom_point(data=infect_predict,aes(x=site,y=estimate__,group=infection),position = position_dodge(width=0.7),shape=23,size=5,fill="white") + scale_fill_manual(values=sitecols[c(3,4)]) + guides(fill="none") + guides(shape = guide_legend(override.aes = list(color = "black") ) )+ theme(legend.position = "right")
rich_infection3 
ggsave2('Richness by Infection Dec23.pdf',rich_infection3,height=15,width=18,units="cm")pH scatterplot
newt1_rich_meta<-mutate(newt1_rich_meta,infection=factor(infection))
richness_ph_scatter<-ggplot(data=newt1_rich_meta, aes(x=as.numeric(as.character(ph)), y=Observed)) +
geom_jitter(aes(fill=site,shape=infection),color="white", size=5.5, width = 0.05) +
theme(text = element_text(size = 40)) +
theme_bw(base_size=15) +
xlab("pH") +
ylab("Alpha Diversity \n (Richness)") + scale_shape_manual(name="Amphibiothecum \n Infection",values=c(21,24),labels=c("No Visible \n Infection","Infected")) + scale_fill_manual(values=sitecols) + labs(fill="Site")
richness_ph_scatter2<- richness_ph_scatter + theme(legend.position = "right") + scale_x_continuous(n.breaks=8)
richness_ph_scatter3<- richness_ph_scatter2 + geom_ribbon(data=ph_predict,aes(x=ph,ymin=lower__,ymax=upper__),alpha=0.3) + geom_line(data=ph_predict,aes(x=ph,y=estimate__),col="blue") + guides(shape = guide_legend(override.aes = list(color = "black"))) + guides(fill = guide_legend(override.aes = list(shape=21))) + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
richness_ph_scatter3
ggsave2('Richness ph Dec23.pdf',richness_ph_scatter3,height=12,width=20,units="cm")Combined Alpha Diversity Plot
#Column Format
alpha_column1<-plot_grid(richness_ph_scatter3,rich_infection3,ncol=1,labels="AUTO",label_size = 30)
ggsave2('Fig_3_Newt Richness Combined Plot.pdf',alpha_column1,width=22,height=40,units="cm")
ggsave2('Fig_3_Newt Richness Combined Plot.tiff',alpha_column1,width=22,height=30,units="cm")BETA DIVERSITY
##Vegan function to convert phyloseq obj. into something Vegan can call directly
vegan_otu <- function(physeq) {
OTU <- otu_table(physeq)
if (taxa_are_rows(OTU)) {
OTU <- t(OTU)
}
return(as(OTU, "matrix"))
}
####TIME POINT 1 PCA####
##CLR Transform - Round 1 sampling only
newt1_clr <- microbiome::transform(newt1, "clr")
#head(otu_table(newt1_clr))
##Extract Matrix and Sample Data
newt1_clr_v<-vegan_otu(newt1_clr) #???
newt1_clr_s<-as(sample_data(newt1_clr),"data.frame")
##Merge C/D infection codes
newt1_clr_s$infectioncode[is.na(newt1_clr_s$infectioncode)] <- "0"
newt1_clr_s$infectioncode<-as.factor(newt1_clr_s$infectioncode)
levels(newt1_clr_s$infectioncode)<-c("0","A","B","C/D","C/D")
levels(newt1_clr_s$infectioncode)## [1] "0" "A" "B" "C/D"###Principal Components Analysis
newt1_pc<-prcomp(newt1_clr_v)
#biplot(newt1_pc)
#plot(newt1r_pc)
summary(newt1_pc)$importance[,1:2]## PC1 PC2
## Standard deviation 15.65805 14.27488
## Proportion of Variance 0.18631 0.15485
## Cumulative Proportion 0.18631 0.34116 newt1_pc_scores<-scores(newt1_pc)
newt1_pc_scores_sub<-newt1_pc_scores[,1:2]
newt1_pc_scores_sub<-cbind(newt1_pc_scores_sub,newt1_clr_s)
##Housekeeping + Label Setup
newt1_pc_scores_sub$site<-as.factor(newt1_pc_scores_sub$site)
newt1_pc_scores_sub$infectioncode<-as.factor(newt1_pc_scores_sub$infectioncode)
### Plot
newt1plotdata<-gg_ordiplot(newt1_pc,groups=newt1_clr_s$site,spider=T,plot=T)
newt1plotdata2<-newt1plotdata$df_ord
newt1plotdata3<-cbind(newt1plotdata2,newt1_clr_s)
newt1plotdata3$infection<-factor(newt1plotdata3$infection)
gg1<-ggplot(newt1plotdata3,aes(x=x,y=y)) + geom_segment(data = newt1plotdata$df_spiders, aes(x = cntr.x, xend = x, y = cntr.y, yend = y, color = Group),
show.legend = FALSE) + geom_point(aes(shape=infection,fill=site),size=7,color="white") + scale_shape_manual(values=c(21,24),name="Amphibiothecum Infection",labels=c("No Visible Infection","Infected"))
gg2<- gg1 + geom_path(data = newt1plotdata$df_ellipse, aes(x = x, y = y, color = Group), show.legend = FALSE)
gg3<- gg2 + theme_bw(base_size=15) + theme(legend.position = "right")
gg4<- gg3 + guides(fill=guide_legend(nrow=3,byrow=TRUE),shape = guide_legend(override.aes = list(color = "black") ),fill= guide_legend(override.aes = list(color = "black",shape=21))) + labs(x="PC1 (18.6%)",y="PC2 (15.8%)",fill="Site") + guides(fill = guide_legend(override.aes = list(shape=21))) ## Warning: Duplicated aesthetics after name standardisation: fill gg5 <- gg4 + scale_fill_manual(values=sitecols) + scale_color_manual(values=sitecols)
gg5
ggsave2('Newt Visit1 CLR PCA.pdf',gg5,width=20,height=15,units="cm")PERMANOVA
##Run PERMANOVA on CLR transformed data
raw_permanova<- adonis2(newt1_clr_v ~ site + infection,data=newt1_clr_s,nperm=999,method="euclidean",by="mar")
raw_permanova #site: r2=26.6% , infection: r2=2.8% ## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = newt1_clr_v ~ site + infection, data = newt1_clr_s, method = "euclidean", by = "mar", nperm = 999)
## Df SumOfSqs R2 F Pr(>F)
## site 3 42012 0.26605 16.0881 0.001 ***
## infection 1 4449 0.02817 5.1112 0.001 ***
## Residual 116 100974 0.63943
## Total 120 157914 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Stacked Barplots
##Filter Samples With No Data using phyloseq object containing clean dataset, which has been filtered from negative/positive control + contaminants, but not filtered for low freq reads.
# ps_clean_filter_filter<-prune_samples(sample_sums(ps_clean_filter)>0,ps_clean_filter)
# newt_phylum <- ps_clean_filter_filter %>% aggregate_taxa(level = "Phylum") %>% microbiome::transform(transform = "compositional")
#Aggregate To Order Level and Transform to Relative Abundance (fraction of reads, rather than count of reads)
#newt_order <- newt.clean.bacteria.sub %>% aggregate_taxa(level = "Order") %>% microbiome::transform(transform = "compositional")
#Phylum Colours
##Plot Again But With Top 'N' taxa
newt_phylumcollapse<- newt1 %>% aggregate_taxa(level="Phylum")
newt_top5phyla = names(sort(taxa_sums(newt_phylumcollapse), TRUE)[1:5]) #What Are the Names of the most abundant phyla?
newt_top5phylum_filter<-subset_taxa(newt1,Phylum %in% newt_top5phyla) #Subset the phyloseq object to those phyla
infectionlong<-ifelse(sample_data(newt_top5phylum_filter)$infection==0,"Uninfected","Infected") #Site and Infection status- Phylum
sample_data(newt_top5phylum_filter)$site_infection<-with(sample_data(newt_top5phylum_filter),paste(site,infectionlong,sep=" "))
##Remake our graph but with grouping by site
newt_top5phylum_plot <- newt_top5phylum_filter %>%
aggregate_taxa(level = "Phylum") %>%
microbiome::transform(transform = "compositional") %>%
plot_composition(sample.sort = "Firmicutes", average_by = "site_infection",group_by = "site")
newt_top5phylum_plot <- newt_top5phylum_plot + theme_bw() +
scale_fill_manual("Phylum", values=moma.colors("Liu",5)) +
theme(legend.position = "bottom") +
theme(text = element_text(size = 22)) + labs(x="Site & Infection Status") +
theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1)) + guides(fill=guide_legend(nrow=3,byrow=TRUE))
newt_top5phylum_plot 
ggsave2('Time 1 Phylum Barplot.pdf',newt_top5phylum_plot,height=20,width=20,units="cm")CCA Model
Model
######
#Extract Matrix and Sample Data
newt1_rare_v<-vegan_otu(newt1)
newt1_rare_s<-as(sample_data(newt1),"data.frame")
##### INDIVIDUAL PREDICTORS
newt_cca<-cca(newt1_rare_v ~ infection + ph + svl,data=newt1_rare_s,strata=frog_rare_s$frog_id)
car::Anova(newt_cca)## Warning in Anova.default(newt_cca): there are rows/columns in vcov.
## that are not in the model matrix:
## CCA1:infection, CCA1:ph, CCA1:svl, CCA2:infection, CCA2:ph, CCA2:svl, CCA3:infection, CCA3:ph, CCA3:svl
## tests may be incorrect## Analysis of Deviance Table (Type II tests)
##
## Response: newt1_rare_v
## Df Chisq Pr(>Chisq)
## infection 1 77.0067 <2e-16 ***
## ph 1 446.2242 <2e-16 ***
## svl 1 0.0228 0.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CCA PLot
#Extract Biplot Dat aand Filter Out the Social Group
bip <- data.frame(scores(newt_cca, choices = 1:2, display = "bp"))
bip_filter<- bip #(bip[-grep("site",rownames(bip)),])
#Scaling Factor
(mul <- ordiArrowMul(bip,display="sites", fill = 0.8))## [1] 0.8298887 bip.scl <- bip_filter * 2
labs<-c("Infection","pH","SVL")
#Extract Plotting Data
newt_cca_plotdata<- ggvegan::fortify(newt_cca) %>% filter(score=="sites")
#Copy Across Metadata
newt_cca_plotdata$site_id<-newt1_rare_s$site
newt_cca_plotdata$infection<-ifelse(newt1_rare_s$infection==0,"No Visible Infection","Infected")
#Plot
newt_cca_plot1<- ggplot(newt_cca_plotdata,aes(x=CCA1,y=CCA2)) + geom_point(color="white",aes(shape=infection,fill=site_id),size=7) + theme_bw(base_size=15) + scale_shape_manual(values=c(24,21)) + guides(fill = guide_legend(override.aes = list(shape = 21))) + scale_fill_manual(values=sitecols) + scale_color_manual(values=sitecols)
#Add Biplot Arrows
newt_cca_plot2<- newt_cca_plot1 + geom_segment(data=bip.scl,aes(x=0,y=0,xend=CCA1,yend=CCA2),arrow=arrow(length=unit(0.2,"cm")),color="grey") #+ lims(x=c(-6,5))
newt_cca_plot3 <- newt_cca_plot2 + geom_text_repel(data=bip.scl,aes(x=CCA1,y=CCA2,label=labs),color="grey40",size=5) + labs(fill="Site",shape="Amphibiothecum Infection") + guides(shape = guide_legend(override.aes = list(color = "black") ))
newt_cca_plot3
ggsave2('Newt Visit1 CCA.pdf',newt_cca_plot3,width=22,height=15,units="cm")Combined Plot
betaplot1<-plot_grid(gg5,newt_cca_plot3,nrow=2,labels="AUTO",rel_widths=c(1,1),label_size = 30)
betaplot1
ggsave2('Newt Visit1 CLR CCA Combined.pdf',betaplot1,width=25,height=25,units="cm")Core Microbiome Visit 1
Newt Site Specific Core
#Unique Sites
newtsites <- unique(as.character(meta(newt1)$site))
#Taxonomy Dataframe
newt1_taxdata<-as.data.frame(tax_table(newt1))
#Write a for loop to go through each of the disease_states one by one and combine identified core taxa into a list.
sites_core <- c() # an empty object to store information
sites_core_taxonomy<-c()
for (n in newtsites){ # for each variable n in DiseaseState
#print(paste0("Identifying Core Taxa for ", n))
ps.sub <- subset_samples(newt1, site == n) # Choose sample from n
core_m <- core_members(ps.sub, # ps.sub is phyloseq selected with only samples from g
prevalence = 0.9)
print(paste0("No. of core taxa in ", n, " : ", length(core_m))) # print core taxa identified in each DiseaseState.
sites_core[[n]] <- core_m # add to a list core taxa for each group.
sites_core_taxonomy[[n]]<-subset(newt1_taxdata,rownames(newt1_taxdata) %in% core_m)
#print(list_core)
}## [1] "No. of core taxa in DP : 27"
## [1] "No. of core taxa in R10 : 10"
## [1] "No. of core taxa in R4 : 6"
## [1] "No. of core taxa in Rnew : 3"#Check Core Taxa Not a FUnction of Sample size
site_ncore_taxa<-data.frame(ntaxa=sapply(sites_core,length))
site_nsamples<-data.frame(nsample=table(meta(newt1)$site))
site_nsamples$ncore<-site_ncore_taxa[,1][match(site_nsamples[,1],rownames(site_ncore_taxa))]
ggplot(site_nsamples,aes(x=nsample.Freq,y=ncore)) + geom_smooth(method = "lm") + geom_point(shape=21,size=5,fill="white") + theme_bw()## `geom_smooth()` using formula = 'y ~ x'
site_ncore_taxa$site<-rownames(site_ncore_taxa)
### Strip Out Taxonomy
site_core_family_prop<-lapply(lapply(sites_core_taxonomy,function(x) table(x$Family)),function(x) x/sum(x))
site_core_family_prop2<-lapply(site_core_family_prop,function(x) data.frame(family=names(x),prop=x))
for(i in 1:length(site_core_family_prop2)){site_core_family_prop2[[i]]$site<-names(site_core_family_prop2)[i]}
site_core_df<-do.call(rbind.data.frame, site_core_family_prop2)
colnames(site_core_df)<-c("Family","Prop","Proportion","Site")
corefamtab<-(table(site_core_df$Family)); corefams<-names(corefamtab)[corefamtab>1]
site_core_df$Family_collapse<-ifelse(site_core_df$Family %in% corefams,site_core_df$Family,"Other")
site_core_df$Family_collapse<-factor(site_core_df$Family_collapse,levels=c("Burkholderiaceae","Comamonadaceae","Oxalobacteraceae","Pseudoalteromonadaceae","Other"))
### Plot
site_core_cols<-c(brewer.pal(5,"Paired"))
site_core_plot1<-ggplot(site_core_df) + geom_bar(aes(x=Site,y=Proportion,fill=Family_collapse),stat="identity",color="grey40") + scale_fill_manual(values=site_core_cols) + theme_bw(base_size=15) + labs(fill="Family") + theme(legend.position = "top") #+ guides(fill=guide_legend(nrow=2))
site_core_plot2<-site_core_plot1 + geom_text(data=site_ncore_taxa,aes(x=site,y=1,label=ntaxa),vjust=-0.2,cex=8) + guides(fill=guide_legend(nrow=3))
site_core_plot2
ggsave2('Core Visit 1 by Site.pdf',site_core_plot2,width=30,height=18,units="cm")Disease Venn
### Plot Venn
core_venn_cols<-brewer.pal(5,"Set2")[c(1:4)]
core_vennplot1<-ggvenn(sites_core,fill_color = core_venn_cols,text_size = 8,set_name_size = 10)
core_vennplot1
Combined Core Plot
core_grid<-plot_grid(site_core_plot2,core_vennplot1,labels="AUTO",ncol=2,label_size = 30)
ggsave2('Fig_6_Newt Core Microbiomes.pdf',core_grid,width=45,height=25,units="cm")Genus Level Core Taxonomy
### Strip Out Taxonomy
site_core_genus<-lapply(lapply(sites_core_taxonomy,function(x) table(x$Genus)),function(x) x/sum(x))
site_core_genus2<-lapply(site_core_genus,function(x) data.frame(genus=names(x),prop=x))
for(i in 1:length(site_core_genus2)){site_core_genus2[[i]]$site<-names(site_core_genus2)[i]}
site_core_genus_df<-do.call(rbind.data.frame, site_core_genus2)
colnames(site_core_genus_df)<-c("Genus","Genus","Proportion","Site")
write.csv(site_core_genus_df,'Site Core ASVs Genus.csv')GENERAL LINEAR LATENT VARIABLE MODELS
Models
#Infected Sites
newt_infected<-prune_samples(sample_data(newt1)$site %in% c("R4","Rnew"),newt1)
## OTUs Top 50
newt_infected_genus<-aggregate_top_taxa2(newt_infected, "Genus", top = 50)
clr_scaled <-microbiome::transform(newt_infected_genus, transform = "clr")
#Extract
ys <- data.frame(t(otu_table(clr_scaled)))
names(ys) <-taxa_names(clr_scaled)
#Predictors
Xs<-data.frame(sample_data(clr_scaled)) %>% dplyr::select(site,infection,svl)
## Model
fit_reduced_scaled <- gllvm(ys, Xs,
num.lv = 2,
formula = ~ site * infection,
family = "gaussian",starting.val='zero')
coefplot(fit_reduced_scaled)
#plot(fit_reduced_scaled)Plots of Effects
#Estimates
df<-coef(fit_reduced_scaled)
est_df<-data.frame(df$Intercept)
est_df2<-data.frame(df$Xcoef)
est_df3<-merge(est_df, est_df2, by = 0)
#Order genera
row.names(est_df3)<-est_df3$Row.names
est_df3<-est_df3[colnames(ys),]
names(est_df3)[1]<- "Genus"
names(est_df3)[2]<- "Intercept"
### COnfidence Intervals
confint_df<-data.frame(confint(fit_reduced_scaled))
#Identify Rows with Main Effects
one_comma<-sapply(rownames(confint_df),function(x) length(gregexpr(":", x, fixed = TRUE)[[1]]))==1
###Strip Out Individual Datasets
##R4
r4_main<-grepl("^Intercept",rownames(confint_df))
int_data_r4<-cbind(est_df3[,c("Genus","Intercept")],confint_df[r4_main,])
##Rnew
rnew_main<-grepl("^Xcoef.siteRnew",rownames(confint_df)) & one_comma==TRUE
int_data_rnew<-cbind(est_df3[,c("Genus","siteRnew")],confint_df[rnew_main,])
##Infection
infection_main<-grepl("Xcoef.infection",rownames(confint_df))
int_data_infection<-cbind(est_df3[,c("Genus","infection")],confint_df[infection_main,])
#Rnew:Infection
rnew_infect_main<-grepl("^Xcoef.siteRnew:infection",rownames(confint_df))
int_data_rnew_infect<-cbind(est_df3[,c("Genus","siteRnew.infection")],confint_df[rnew_infect_main,])
#Extra Variables and Rename Columns
colnames(int_data_rnew)<-c("Genus","Estimate","l95","u95")
colnames(int_data_infection)<-c("Genus","Estimate","l95","u95")
colnames(int_data_rnew_infect)<-c("Genus","Estimate","l95","u95")
int_data_rnew$trait<-"Site Rnew"
int_data_infection$trait<-"Infection"
int_data_rnew_infect$trait<-"Rnew:Infection"
newt_mod_plotdata<-rbind(int_data_infection,int_data_rnew_infect)
#Order
newt_mod_plotdata$trait<-factor(newt_mod_plotdata$trait,levels=c("Infection","Rnew:Infection"))
newt_mod_plotdata2<- newt_mod_plotdata %>% group_by(trait) %>% arrange(Estimate,.by_group=T)
newt_mod_plotdata2$Genus<-factor(newt_mod_plotdata2$Genus,levels=unique(newt_mod_plotdata2$Genus))
#Significance
newt_mod_plotdata2$Sig<- !data.table::between(0, newt_mod_plotdata2$l95, newt_mod_plotdata2$u95)
sig_col<-brewer.pal(5,"Set1")[2]
### plot
newt_plot1<-ggplot(newt_mod_plotdata2,aes(x=Estimate,y=Genus)) + geom_errorbar(aes(y=Genus,xmin=l95,xmax=u95,color=Sig),alpha=0.6) + geom_point(size=5,shape=21,color="gray40",aes(fill=Sig),alpha=0.9)
newt_plot2<- newt_plot1 + theme_bw(base_size = 15) + geom_vline(xintercept=0,linetype="dashed") + scale_color_manual(values=c("gray40",sig_col)) + scale_fill_manual(values=c("white",sig_col)) + guides(fill="none",color="none") + facet_wrap(.~trait,scales = "free_x")
newt_plot2
Significant Bacterial Trends
##Pseudomonas
pseudo_df<-data.frame(count=as.numeric(otu_table(clr_scaled)[which(rownames(otu_table(clr_scaled))=="Pseudomonas"),]),site=sample_data(newt_infected_genus)$site,infection=sample_data(newt_infected_genus)$infection)
pseudo_df$infection2<-ifelse(pseudo_df$infection==1,"Infected","No Visible Infection")
pseudo_plot1<-ggplot(pseudo_df,aes(x=site,y=count,group=infection2)) + geom_point(aes(fill=site,shape=infection2),size=5,color="gray77",position = position_jitterdodge(dodge.width = 0.75)) + theme_bw(base_size = 18) + scale_fill_manual(values = sitecols[3:4]) + scale_shape_manual(values=c(24,21)) + guides(fill="none") + labs(x="Site",y="Count (CLR)",shape="Amphibiothecum \nInfection") + theme(legend.position = "top") + stat_summary(fun = "mean",geom = "point",
col = "white",size = 5,shape = 23,fill = "lightblue",position = position_dodge(width = 0.75)) + geom_text(x=1.5,y=5,label="Pseudomonas",cex=8) + guides(shape=guide_legend(nrow=2,byrow=TRUE))
pseudo_plot1 
##Bacteroides
bacteroides_df<-data.frame(count=as.numeric(otu_table(clr_scaled)[which(rownames(otu_table(clr_scaled))=="Bacteroides"),]),site=sample_data(newt_infected_genus)$site,infection=sample_data(newt_infected_genus)$infection)
bacteroides_df$infection2<-ifelse(bacteroides_df$infection==1,"Infected","No Visible Infection")
bacteroides_plot1<-ggplot(bacteroides_df,aes(x=site,y=count,group=infection2)) + geom_point(aes(fill=site,shape=infection2),size=5,color="gray77",position = position_jitterdodge(dodge.width = 0.75)) + theme_bw(base_size = 18) + scale_fill_manual(values = sitecols[3:4]) + scale_shape_manual(values=c(24,21)) + guides(fill="none") + labs(x="Site",y="Count (CLR)",shape="Amphibiothecum \nInfection") + theme(legend.position = "top") + stat_summary(fun = "mean",geom = "point",
col = "white",size = 5,shape = 23,fill = "lightblue",position = position_dodge(width = 0.75)) + guides(shape=guide_legend(nrow=2,byrow=TRUE)) + geom_text(x=1.5,y=4,label="Bacteroides",cex=8)
bacteroides_plot1 
##Janinthobacterium
janinthobacterium_df<-data.frame(count=as.numeric(otu_table(clr_scaled)[which(rownames(otu_table(clr_scaled))=="Janthinobacterium"),]),site=sample_data(newt_infected_genus)$site,infection=sample_data(newt_infected_genus)$infection)
janinthobacterium_df$infection2<-ifelse(janinthobacterium_df$infection==1,"Infected","No Visible Infection")
janinthobacterium_plot1<-ggplot(janinthobacterium_df,aes(x=site,y=count,group=infection2)) + geom_point(aes(fill=site,shape=infection2),size=5,color="gray77",position = position_jitterdodge(dodge.width = 0.75)) + theme_bw(base_size = 18) + scale_fill_manual(values = sitecols[3:4]) + scale_shape_manual(values=c(24,21)) + guides(fill="none") + labs(x="Site",y="Count (CLR)",shape="Amphibiothecum \nInfection") + theme(legend.position = "top") + stat_summary(fun = "mean",geom = "point",
col = "white",size = 5,shape = 23,fill = "lightblue",position = position_dodge(width = 0.75)) + guides(shape=guide_legend(nrow=2,byrow=TRUE)) + geom_text(x=1.5,y=3,label="Janthinobacterium",cex=8)
janinthobacterium_plot1 
##Chryseobacterium
chryseobacterium_df<-data.frame(count=as.numeric(otu_table(clr_scaled)[which(rownames(otu_table(clr_scaled))=="Chryseobacterium"),]),site=sample_data(newt_infected_genus)$site,infection=sample_data(newt_infected_genus)$infection)
chryseobacterium_df$infection2<-ifelse(chryseobacterium_df$infection==1,"Infected","No Visible Infection")
chryseobacterium_plot1<-ggplot(chryseobacterium_df,aes(x=site,y=count,group=infection2)) + geom_point(aes(fill=site,shape=infection2),size=5,color="gray77",position = position_jitterdodge(dodge.width = 0.75)) + theme_bw(base_size = 18) + scale_fill_manual(values = sitecols[3:4]) + scale_shape_manual(values=c(24,21)) + guides(fill="none") + labs(x="Site",y="Count (CLR)",shape="Amphibiothecum \nInfection") + theme(legend.position = "top") + stat_summary(fun = "mean",geom = "point",
col = "white",size = 5,shape = 23,fill = "lightblue",position = position_dodge(width = 0.75)) + guides(shape=guide_legend(nrow=2,byrow=TRUE)) + geom_text(x=1.5,y=0,label="Chryseobacterium",cex=8)
chryseobacterium_plot1 
Combined Model and Output Plot
#Abundance Grid
gllvm_abund_grid<-plot_grid(chryseobacterium_plot1,janinthobacterium_plot1,bacteroides_plot1,pseudo_plot1,nrow=2,
label_size = 20)
#Stitch onto Model
gllvm_combined<-plot_grid(newt_plot2,gllvm_abund_grid,nrow=2,labels="AUTO",label_size=30)
ggsave2('Newt GLLVM Combined Output.png',gllvm_combined,width=30,height=50,units="cm")Correlation Plots
#### Correlations
cr1<-data.frame(getResidualCor(fit_reduced_scaled))#
names(cr1)<-names(ys)
names(cr1)<-abbreviate(names(cr1), minlength = 15)
rownames(cr1)<-abbreviate(rownames(cr1), minlength = 15 )
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:stats':
##
## filter #devtools::install_github("kassambara/ggcorrplot")
library(ggcorrplot)##
## Attaching package: 'ggcorrplot'## The following object is masked from 'package:rstatix':
##
## cor_pmat #install.packages("ggpubr")
library(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legend cr2<-cor_pmat(cr1)
corplot<-ggcorrplot(cr1,
hc.order = TRUE,
outline.col = "white",
type = "full",
ggtheme = ggplot2::theme_minimal(base_size = 10),
tl.cex = 12,
p.mat = cr2,
sig.level = 0.05,
lab_size = 15,
#show.diag = F,
insig = "blank",
# colors = c("#6D9EC1", "white", "#E46726"))
colors = c("blue", "white", "red"))+
theme(axis.text.x = element_text(angle = 90, vjust=0.5))+
theme(plot.margin=unit(c(0.2,0.2,0.2,2),"cm"))## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggcorrplot package.
## Please report the issue at <https://github.com/kassambara/ggcorrplot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated. corplot
## Save Grid
#newt_mod_combined1<-ggarrange(newt_plot2,corplot,nrow=2,widths=c(1,0.8),align="h",labels="AUTO",font.label=list(size=30))
#newt_mod_combined1
cowplot::ggsave2('Newt GLLVM Residual Correlations.pdf',corplot,width=25,height=25,units = "cm")PREDICTED FUNCTIONAL ANALYSIS
Prepare for MicFunPred
Adapted from https://forum.qiime2.org/t/importing-dada2-and-phyloseq-objects-to-qiime-2/4683
# ## Trim Newt1
# newt1_filter<-prune_taxa(taxa_sums(newt1)>200,newt1)
#
# #Make Fasta File From Sequences
# #remotes::install_github("alexpiper/seqateurs")
# library(seqateurs)
# ps_to_fasta(newt1_filter)
#
# #Export Feature Table
# #library(biomformat);packageVersion("biomformat")
# otu<-t(as(otu_table(newt1_filter),"matrix")) # 't' to transform if taxa_are_rows=FALSE
# #Rename ASVs to Match FASTA above
# rownames(otu)<-paste0("ASV_",seq(nrow(otu)))
# #Fix Sample Names
# true_samples_newt1filter<-colnames(otu)
# colnames(otu)<-paste0("sample",seq(ncol(otu)))
# #Write File
# # otu_biom<-make_biom(data=otu)
# # write_biom(otu_biom,"Newt1_filter.biom")
# write.table(otu,'newt1_filter.tsv',sep="\t",row.names=TRUE, col.names=TRUE, quote=FALSE)
#
# ## Sample Data
# newt1_samplemeta<-as(sample_data(newt1_filter),"data.frame")
# write.csv(newt1_samplemeta,"Newt1_Filter_sample-metadata.csv")BETA DIVERSITY OVER TIME
# Visit 1 and 2
newt2<-prune_samples(sample_data(ps_clean_filter)$visit %in% seq(2) & !sample_data(ps_clean_filter)$site %in% c("R44","Rnew"),ps_clean_filter)
##CLR Transform - Round 1 & 2 sampling only
newt2_clr <- microbiome::transform(newt2, "clr")
##Extract Matrix and Sample Data
newt2_clr_v<-vegan_otu(newt2_clr) #???
newt2_clr_s<-as(sample_data(newt2_clr),"data.frame")
###Principal Components Analysis
newt2_pc<-prcomp(newt2_clr_v)
summary(newt2_pc)$importance[,1:2]## PC1 PC2
## Standard deviation 15.37126 14.99021
## Proportion of Variance 0.16088 0.15301
## Cumulative Proportion 0.16088 0.31389## Metadata
newt2_pc_scores<-scores(newt2_pc)
newt2_pc_scores_sub<-newt2_pc_scores[,1:2]
newt2_pc_scores_sub<-cbind(newt2_pc_scores_sub,newt2_clr_s)
##Housekeeping + Label Setup
newt2_pc_scores_sub$site<-as.factor(newt2_pc_scores_sub$site)
newt2_pc_scores_sub$visit<-as.factor(newt2_pc_scores_sub$visit)
newt2_pc_scores_sub$site_visit<-with(newt2_pc_scores_sub,paste(site,visit))
### Plot Data
newt2plotdata<-gg_ordiplot(newt2_pc,groups=newt2_pc_scores_sub$site_visit,spider=T,plot=T)
#newt2plotdata$plot + scale_color_brewer(palette = "Paired") + theme_bw() + geom_point(size=5)
newt2plotdata2<-newt2plotdata$df_ord
newt2plotdata3<-cbind(newt2plotdata2,newt2_clr_s)
newt2plotdata3$Group<-factor(newt2plotdata3$Group)
newt2plotdata3$infection_text<-ifelse(newt2plotdata3$infection==0,"No Visible Infection","Infected")
gg_twovisits1<-ggplot() + geom_segment(data = newt2plotdata$df_spiders, aes(x = cntr.x, xend = x, y = cntr.y, yend = y, color = Group),
show.legend = FALSE) + geom_point(data=newt2plotdata3,aes(x=x,y=y,fill=Group,shape=infection_text),color="white",size=5)
gg_twovisits2<- gg_twovisits1 + geom_path(data = newt2plotdata$df_ellipse, aes(x = x, y = y, color = Group), show.legend = FALSE)
gg_twovisits3<- gg_twovisits2 + theme_bw(base_size=15) + theme(legend.position = "right") + scale_shape_manual(values=c(24,21))
gg_twovisits4<- gg_twovisits3 + guides(shape = guide_legend(override.aes = list(color = "black") ) ) + labs(x="PC1 (16.09%)",y="PC2 (15.3%)",fill="Site & Visit",shape="Amphibiothecum \n Infection") + guides(fill = guide_legend(override.aes = list(shape=21)))
gg_twovisits5<-gg_twovisits4 + scale_color_brewer(palette="Paired") + scale_fill_brewer(palette = "Paired")
gg_twovisits5
ggsave2('Newt 2 Time Point PCa.pdf',gg_twovisits5,width=24,height=15,units="cm")PERMANOVA
##Run PERMANOVA on CLR transformed data
visit_permanova<- adonis2(newt2_clr_v ~ site*visit + infection,data=newt2_clr_s,nperm=999,method="euclidean")
visit_permanova ### Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = newt2_clr_v ~ site * visit + infection, data = newt2_clr_s, method = "euclidean", nperm = 999)
## Df SumOfSqs R2 F Pr(>F)
## Model 6 61508 0.37731 10.604 0.001 ***
## Residual 105 101509 0.62269
## Total 111 163017 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1